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Two fermions occupying the same site of a lattice model with strongly repulsive Hubbard-type 
interaction form a doublon, a long-living excitation the decay of which is suppressed because of 
energy conservation. By means of an exact-diagonalization approach based on the Krylov-space 
technique, we study the dynamics of a single doublon, of two doublons, and of a doublon in the 
presence of two additional fermions prepared locally in the initial state of the extended Hubbard 
model. The time dependence of the expectation value of the double occupancy at the different 



fS| sites of a large one-dimensional lattice is analyzed by perturbative arguments. In this way the 

spatiotemporal evolution of the doublon, its decay on a short time scale and the long-time average 
of the total double occupancy can be understood. We demonstrate how the dynamics of a doublon in 
the initial state is related to the spectrum of two-fermion excitations obtained from linear-response 
^ theory, we work out the difference between doublons composed of fermions vs. doublons composed of 

^ bosons and show that despite the increase of phase space for inelastic decay processes, the stability 

of a doublon is enhanced by the presence of additional fermions on an intermediate time scale. 

CO PACS numbers: 71.10.Fd, 67.85.-d, 78.47.D- 
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I. INTRODUCTION 

Since the seminal work of Jaksch et al.,^ ultracold 
atomic gases in optical lattices have served as a valuable 
testing ground for the rich phenomenology of many-body 
models which originally were introduced in the context 
of condensed- matter physics. A nice example is the 
concept of repulsively bound pairs of fermions which can 
be studied in the strong-coupling regime of the Hubbard 
model or, as shown recently,^ in an organic salt at room 
temperature by means of ultrafast optical spectroscopy. 
Repulsively bound pairs, named doublons, are already 
known since the early work of Hubbard^ and were lately 
addressed in both theoretical and experimental work in 
bosonic^^^ as well as fermionic^^^^ Hubbard-type mod- 
els. The fermionic case directly refers to condensed- 
matter systems, such as strongly correlated electrons in a 
valence band of transition-metals and their compounds, 
and two-particle electron spectroscopy. 

A doublon is a pair of two fermions tightly bound to 
each other. The pair is itinerant. It propagates through 
the lattice and thereby acquires a certain energy disper- 
sion. The pair may decay into its constituents. However, 
for strongly repulsive interaction [/ > 0, this decay is 
suppressed very efficiently. The stability of the doublon 
appears as counterintuitive since an energy of the order of 

> would gained if the two fermions were propagating 
through the lattice independently. There is, however, a 
"repulsive binding" originating from energy conservation: 
For U much larger than the nearest-neighbor hopping J, 
the excess energy U cannot be accommodated in the ki- 
netic energy of the two independent fermions which at 
most amounts to twice the bare bandwidth W oc J. 

In the strong-coupling limit, doubly occupied sites are 
created in different types of electron spectroscopies:^^ 
The spectral function Aij{uj) at positive frequencies uo > 
0, obtained from the imaginary part of the one-particle 
Green's function ((c-^;c^^)), is related to inverse pho- 



toemission, and the upper Hubbard band in the spec- 
tral function describes a final state with doubly occupied 
sites. The lower Hubbard band represents the analogue 
of the upper Hubbard band in case of photoemission. For 
the Hubbard model on a bipartite lattice at half-filling, it 
is obtained from the upper one by a particle-hole trans- 
formation and thus describes repulsively bound holes. 
Doublons can also created in an otherwise empty va- 
lence band in a two-particle process, such as appear ance- 
potential-spectroscopy (APS), i.e. the "time-inverse" of 
Auger-electron spectroscopy (AES). Here, two additional 
electrons (holes) are created, preferably at the same site, 
in the final state of APS (AES). Doublon bound states in 
APS/ AES show up in the local two-particle Green's func- 
tion {{c-^c--; cj^cj-)) as is well known from Cini-Sawatzky 
theory. ^^^^ Furthermore, doublons appear in particle- 
hole excitations associated with Green's functions of the 
type ((4cr^jcr5 ^Icr^^cr))- meutioued cases, a doublon 

would be identified with a long-lived excitation at ener- 
gies of the order of U. 

Since a pair of fermions has bosonic character, the 
exciting question arises whether a macroscopically large 
number of doublons could Bose condensate at sufficiently 
low temperatures and high densities. This has been stud- 
ied theoretically for doublons of bosonic^ and of fermionic 
constituents.^^ The two main questions in this context 
concern the doublon stability and the effective interac- 
tion between doublons: First, a sufficiently long lifetime 
of doublons is required for a possible Bose condensation 
taking place in a metastable state. Recent experiments 
with fermionic atoms in optical lattices^^ in fact give a 
lifetime which increases exponentially with U. Second, 
the physical interactions between the constituents give 
rise to an effective doublon-doublon interaction in an 
effective low-energy theory. A strong repulsive [/, for 
example, leads to an effectively attractive interaction be- 
tween doublons formed by two bosons. It has been shown 
that this inhibits condensation but rather favors phase 
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separation. 

The real-time dynamics of a spatially extended system 
of strongly correlated fermions poses a notoriously com- 
plex many-body problem which is hardly accessible to 
exact analytical or numerical methods. With the present 
paper we study a simplified problem with a drastically 
reduced Hilbert space dimension by focussing on a few 
fermions only. The time evolution of this few-fermion 
quantum system is accessible by a numerically exact 
Krylov-space approach. ^^^^ Specifically, we consider sys- 
tems of two and of four spinful fermions on large one- 
dimensional lattices of typically L = 50 — 100 sites with 
nearest-neighbor hopping J and a strong on-site Hubbard 
repulsion /7, and with a nearest-neighbor interaction V 
in addition. The comparatively simple but numerically 
exact Krylov technique allows to study the short- and 
long-time dynamics of a single doublon, of two doublons 
or of a doublon in the background of two additional un- 
bound fermions. 

Our study tries to shed some light on the following 
issues discussed extensively in the recent literature: 

A doublon prepared at a definite site initially, starts to 
propagate through the lattice. The presence of a nearest- 
neighbor repulsion can dramatically change the character 
of this propagation. We show how the resulting propaga- 
tion patterns can be understood in detail with the help of 
an effective low-energy model or perturbative arguments 
in the strong-coupling limit and argue that for certain 
resonance conditions a much higher mobility of the dou- 
blon is possible. For U = V m. particular, the propaga- 
tion is better understood in terms of an extended object. 

The above-mentioned energy-conservation argument is 
somewhat imprecise as even a single doublon in an other- 
wise empty band may undergo a decay into its fermionic 
constituents if U is finite. We will analyze the manifesta- 
tion of "energy conservation" by studying the short-time 
dynamics of a doublon. 

The real-time dynamics of a quantum system in a 
highly excited state on the one hand and the spectrum 
of excitations out of thermal equilibrium, as obtained in 
linear-response theory, on the other hand are usually two 
completely different issues. Here, we discuss a one-to-one 
relation that is obtained for the case of a single doublon. 
This relation and also the method of canonical transfor- 
mation helps to understand the physics of the long-time 
stability of a single doublon, i.e. the U dependence of the 
time-averaged double occupancy. 

The real-time dynamics of two doublons, initially pre- 
pared at a distance is much more complicated. We 
discuss the main trends as functions of U and V with 
the help of perturbative arguments in the strong-coupling 
limit. Particularly the V dependencies are interesting as 
there is a reduced effective doublon-doublon attraction in 
the Fermi opposed to the Bose case which is important 
to understand the competition between Bose condensa- 
tion and phase separation of doublons.^' Our study sug- 
gests that there is no clustering of doublons consisting of 
fermions unless an additional V interaction is present. 



For lattice fermion models with a finite particle den- 
sity, i.e. for a thermodynamically relevant number of 
fermions, the problem becomes intractable by exact nu- 
merical means. One generally expects, however, that 
with the presence of many additional degrees of freedom 
there is an enhanced probability for doublon decay as the 
doublon energy can by accommodated among different 
particles in a high-order scattering event. This is already 
seen by means of the ladder approximation applicable 
to the low-density limit where a strong initial decay at 
short times is observed followed by a slow exponential de- 
cay at long times. For the presently studied case of four 
fermions, one would not expect an exponential decay law 
but a decreased long-time stability of a doublon if pre- 
pared on top a two-fermion ground state. We will show 
and explain, however, that the presence of the additional 
fermion rather leads to an enhanced stability. This find- 
ing is discussed in the context of recent time-dependent 
density-matrix renormalization-group studies. ^^'^^ 

The paper is organized as follows: The next Sec. H 
introduces the model and the Krylov time-evolution ap- 
proach. We start with the analysis of single-doublon 
propagation in Sec. HI, discuss the effects of the nearest- 
neighbor interaction in Sec. IV and the short-time decay 
in Sec. V. The relation to APS is worked out in Sec. VI, 
and the long-time stability is discussed in Sec. VII. The 
second part of the paper is devoted to our four-fermion 
results: We discuss the dynamics of two doublons in Sec. 
VIII and doublon-fermion scattering in Sec. IX. Final re- 
marks and conclusions are given in Sec. X. 

II. KRYLOV APPROACH TO THE EXTENDED 
HUBBARD MODEL 

Ultracold atoms, loaded into an optical lattice are sub- 
ject to different kinds of interaction.^^ In the simplest 
cases these are short-ranged, like van der Waals forces, 
scaling as and hence approximately act on-site only. 
Depending on the atomic species, however, more gen- 
eral interactions can occur. For example, polarized dipo- 
lar atoms experience a dipole-dipole interaction given by 
/7dd oc (1 — 3cos^^)/r^. Depending on the angle 9 be- 
tween the dipole moments and their relative displace- 
ment vector, this can either be repulsive or attractive. 
It is comparatively long-ranged and usually modeled as 
an interaction between nearest neighbors. Overall, this 
motivates the extended Hubbard model: 

^ = --^ E E 4s-. + E ^u^it 

<ij> (7 i 

+ E E^-^5-'=-^-^+^c^+^^' (1) 

<ij> era' 

which also applies as a model description to elec- 
trons interacting via the screened Coulomb repulsion 
in condensed-matter systems, e.g. transition-metal com- 
pounds, if orbital degrees of freedom can be neglected. 
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Here, i and j refer to the sites of a one- or higher- 
dimensional lattice, < ij > denotes nearest neighbors, 
and a =t,i is the spin projection. J is the nearest- 
neighbor hopping, and U and V the on-site and the 
nearest-neighbor interaction strength. To allow for com- 
parison with the results of previous calculations, our cal- 
culations are performed for a one-dimensional lattice. 

Our central object of interest is the time-dependent 
expectation value of both, the local and the total double 
occupancy, namely {T>i{t)) and {V(t)) = {Vi{t)), re- 
spectively. Here, the local double-occupation operator is 
given by Vi = fi^^^n^^ where n^^ = c\^c-^ is the number 

operator and where denotes the annihilation (cre- 
ation) operator for a fermion at i with spin <j. The time 
dependence of the expectation value is due to the time de- 
pendence of the system's state \i^(t)) = ex.p{—il-Lt) |V^ini) 
where |V^ini) is the state in which the system was prepared 
initially at time t = 0. 

For systems with moderately large Hilbert-space di- 
mensions (i, the numerically exact time evolution of a 
given initial state is accessible by means of a time- 
dependent Krylov-space technique. ^^^^ For a given vec- 
tor u the n-th Krylov subspace of the full Hilbert space 
is defined by^^ 

^ri{u, n) := span [u, Uu, . . . , W'^u] . (2) 

Typically, the Krylov-space dimension n <^ d. An or- 
thogonal basis of can be obtained efficiently via the 
Lanczos recursion formula^^ 

Uk^i = Uuk - akUk - h\uk-i (/c = 0, . . . , n - 1) , (3) 

with the coefficients ak = {uk\Huk) / {uk\uk) and h\ = 
{uk\uk) / {uk-i\uk-i) and the initial values bo = and 
U-i = 0. In the normalized Lanczos basis {vi}, with 
Vi = Ui/\\ui\\, the Hamiltonian is represented by a tridi- 
agonal matrix T with diagonal elements ao, . . . , ftn-i and 
secondary diagonal elements 6i, . . . , bn-i- Hence, we can 
write T = V^HV, where the matrix V is made up by the 
basis vectors Vi^ i.e. V = {vq, . . . , Vn-i)- 

The time evolution of a state G = ^n{t) ap- 
proximates its time evolution in the whole Hilbert space: 
ilj{t + At) ^ Ve-^^(^+^^) VtV;(t). Here is chosen to 
be the start vector of the Lanczos recursion (Eq. 3), i.e. 
the Krylov space at time t is adjusted to the system's 
state at t. For a given small time step At, the approxi- 
mation can be controlled to a high accuracy by adjusting 
the Krylov-space dimension. Longer time evolutions are 
carried out successively by using the propagated state 
as the new initial state and adapting T and V after each 
Lanczos time step. It is important to note, that this kind 
of approximation preserves the unitarity of the time evo- 
lution. 

Since the diagonalization of the fairly small n x n ma- 
trix T is numerically cheap, the computational effort 
is dominated by the n — 1 matrix- vector-multiplications 
that are necessary to construct the Lanczos basis and by 
the number of time steps. In this work we dealt with 



Hilbert spaces with d = 10^ . . . 10^ dimensions. For cal- 
culations where e.g. 200 time steps At = 0.5 are per- 
formed, highly accurate results are obtained using Krylov 
spaces with less then n ~ 20 dimensions only. 

In the following we concentrate on a one-dimensional 
lattice with L sites and two or four fermions with equal 
number of up and down spins. Thereby, different pro- 
cesses, such as the propagation and decay of a single 
doublon as well as doublon-fermion and doublon-doublon 
scattering can be studied. For two fermions, the Hilbert- 
space dimension is d = and we opt for a lattice with 
L = 100 sites. For four fermions, it is d = — 1)^/4 

and we shorten the lattice to 50 sites. In either case, 
periodic boundary conditions are assumed. 



III. PROPAGATION OF A SINGLE DOUBLON 

To begin with, we consider the two-fermion system and 
assume that initially, at time t = 0, both fermions are at 
the same site io^ i.e. |V^ini) = cj^^cj^^ |0). Fig. 1 (left part) 
shows the time evolution of the expectation value of the 
local double occupancy at ^ = and for strong on-site 
interaction U = 8 J. The nearest- neighbor hopping J = 1 
fixes the energy and time scales. We notice different ef- 
fects. First of all, the doublon delocalizes. The double 
occupancy {Vi{t)) at the site where the doublon has been 
prepared initially (zq = 50) quickly decreases, and in the 
course of time {Vi{t)) basically spreads out over the en- 
tire lattice. For the time scale t < 100 shown in the 
figure, the "light cones" do not yet interfere through the 
periodic boundary. Second, there is doublon decay. The 
top panel of Fig. 1 shows the total double occupancy 
{V{t)) = {Vi{t)). There is a significant decay from 
the initial value {V{t)) = 1 to about {V{t)) ^ 0.9 in a 
very short time t < 0.5 (not resolved on the scale of the 
figure), followed by an almost constant trend. The tiny 
fiuctuations around the constant "final" value are simply 
refiecting the fact that the total double occupancy does 
not commute with the Hamiltonian. 

Except for the decay of the doublon, all the details 
of the entire propagation profile are fully captured by a 
simple analytical description in an effective low-energy 
model, see the right panel of Fig. 1. To this end, we 
employ the technique of unitary transformations^^^^ to 
project out the energetically well separated high-energy 
part of the spectrum, thereby generating effective low- 
energy couplings perturbatively, in powers of J/U. The 
main steps are described below (see also Refs. 8, 9 and 
17). The goal is to derive a low-energy Hamiltonian pre- 
serving the total double occupancy. 

First, the hopping term Hj is subdivided into parts 
preserving or changing the total double occupancy of the 
system. Expressing the identity by number operators for 
particles and holes, namely lio- = h^^ + n^^, one may 
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FIG. 1. Temporal evolution of the expectation value (see 
color code) of the double occupancy at sites i = 1, L. Left: 
Numerical results for the one-dimensional Hubbard model 
(V — 0) with L — 100 sites and periodic boundary conditions 
at [/ = 8 J. Initially, at t = 0, the two-fermion state has been 
prepared as a doublon at iq = 50. Top: Time evolution of 
the expectation value of the total double occupancy. Right: 
Corresponding analytical results of the effective model, see 
Eq. (9). Note that most of the color range is used to display 
expectation values less than 0.1, as shown in the color bar 
below. The time t is measured in units of the inverse hopping 

i/j. 

write 

= --^ E E ("?^4s>,^^ + k-A.c,^%) 

<ij> cr 

- E E"i-^l-s-^'^ ~ E E'^i-^l-s-^ 

<ij> cr <ij> cr 

=:H'j+Hj + nj, (4) 

where the double occupancy is raised/ lowered by and 
preserved by Hj, since 

[Hu,Wj] = vUWj ve{0,±}. (5) 

The unitary transformation is performed perturbatively: 

W = e'^Ue-'^ ^n^i[S,n]^- [5, [5, H]] + . . . . (6) 

'Hj can be eliminated by choosing the generator to be 
S = - {j (Hj -nj). Up to order J'^/U, we end up with 
the effective model 

H.s = n°j + nu + ^[n+,Hj] , (?) 

which, besides the total particle number, conserves the 
total double occupancy in addition. We can therefore 
restrict ourselves to a system without any singly occupied 
site. This yields, after some algebra, a Hamiltonian, 

<ij> i <'ij> 

given in terms of doublon degrees of freedom only: d] = 
cj^cj^ and d- describe hard-core bosons with the con- 
straint dl'^ = 0. Furthermore, hf = did- = fil^fil^ is 



the local doublon number. Eq. (7) takes the form of an 
extended Bose Hubbard model featuring an (in case of 
positive U) attractive nearest-neighbor interaction where 
J' = AJ'^ /U is the effective hopping parameter. 

For a system with a single doublon only, the interac- 
tion term can be disregarded, and the resulting free tight- 
binding Hamiltonian is diagonalized by Fourier transfor- 
mation. In the limit L ^ oo, the time dependent local 
double occupancy in the effective model is then found to 
be given by the k-ih. Bessel function of the first kind Jk^ 

{Vf{t)) = Jl,^iJ't), (9) 

if the doublon was prepared at site initially. Note 
that the total double occupancy is conserved, since 
Er=-oo^fc(2^) = lforall X. 

The time dependence of the expectation value of the 
local double occupancy, as given by Eq. (9), is shown in 
Fig. 1 (right). While effects due to doublon decay are ne- 
glected at this level, doublon-propagation effects should 
be captured qualitatively correct. Comparing with the 
exact numerical result (Fig. 1, left), we note that the 
effective model provides an excellent description of the 
propagation already for U = SJ. 

The effect of varying U can be seen in Fig. 2. The 
panels Fig. 2(c), (h) and (m) give the result of the full 
model for = 0, = 5 J and U 10 J. We note that 
the mobility of the doublon decreases with increasing U 
which, in the effective model, is due to the reduced dou- 
blon hopping ^ 1/U. The interference pattern visible for 
U = 5J in panel Fig. 2(h) is due to the finite system 
size and periodic boundary conditions. Apart from that, 
however, the pattern does not change much qualitatively 
as compared to U = 10 J. This is worth mentioning since 
U = 5 J is well below the critical U (of the order of twice 
the free bandwidth 2W = 8J) at which the two-particle 
excitation spectrum, related to the APS Green's func- 
tion ((c-^c--; ct^ct-)), does change qualitatively since the 
correlation satellite splits off (see Ref. 25, for example). 
This reminds us that there is a clear conceptual difference 
between the two-particle spectrum that refers to excita- 
tions starting from the system's ground state on the one 
hand and the temporal evolution of a highly excited ini- 
tial state on the other. 



IV. EFFECTS OF NEAREST-NEIGHBOR 
INTERACTION 

The remaining panels of Fig. 2 show propagations 
patterns for finite nearest-neighbor interaction V. For 
U = lOJ, see last row in Fig. 2, we find a decreasing 
mobility of the doublon with increasing difference be- 
tween the on-site and the nearest-neighbor interaction 
strengths U — V. Similar to the discussion in the preced- 
ing section, this trend is easily explained in an effective 
model that preserves the total double occupancy. This 
can be derived, for example, by standard second-order 
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(a) t/ = 0,y = -10 {h)U = 0,V = -5 {c)U = 0,V = {d)U = 0,V = 5 {e)U = 0,V = 10 
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(i)U = 5,V = -10 (g) U = 5,V = -5 {h)U = 5,V = (i)U = 5,V = 5 Q)U = 5,V = 10 
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(k) U = 10,V = -10 (\)U = 10,V = -5 (m)U = 10,V = (n)U = 10,V = 5 (o) U = 10,V = 10 

FIG. 2. Time-dependent expectation value of the local (main panels) and total double occupancy (small top panels) for two 
fermions initially prepared at the same site io of a one-dimensional lattice with L = 100 sites and periodic boundary conditions. 
Results for on-site interaction [/ = 0, 5, 10 and nearest-neighbor coupling V — 0, ±5, ±10, as indicated. Note that for the total 
double occupancy the scale of the ^-axis differs from case to case. The color code is the same as in Fig. 1 and the same for all 
plots. 



perturbation theory around the J = limit and yields 
an effective doublon hopping amplitude 



This corresponds to a sequence of two virtual hopping 
processes: In the first, one of the two fermions composing 
the doublon hops to a nearest neighbor site. Thereby, 
for U > V {U < V^), the energy U — V is gained (has to 
be paid). The second nearest-neighbor hopping process 
leads to the recombination of the doublon, either at the 
same or at one of the adjacent sites. 

Looking at the cone angle of the "light cone" in the 
propagation patterns in the last row and comparing the 
results with - 1/ = 20 to - F = 5, the effective 
description yields the correct trend: As the expectation 
value of the double occupancy in the effective model de- 
pends on the product of J' and t only, see Eq. (9), the 
time axis scales linearly with J^ The effective doublon 
hopping also explains that the patterns in panels Fig. 2(f) 
and (1) and the patterns in (g) and (m) as well as (h) and 
(n) are almost equal as /7 — F is constant respectively. 

For U = 5J, see middle row in Fig. 2, the results for 
U — V = 5 and U — V = —5 differ significantly although 
they should be described by the same effective hopping 
J\ apart from the sign. The sign, however, has no effect. 



The difference is rather due to the residual influence of 
virtual processes of fourth order in J where one of the 
fermions hops two sites away, followed by a recombina- 
tion of the doublon. This leads to an asymmetry between 
the two cases, ±|/7— y|, since for U > V all three interme- 
diate states have lower energy while for U < V two states 
are higher in energy and one lower. With increasing in- 
teraction strengths U and we find this asymmetry to 
be less and less efficient as expected. 

As can be seen by comparing panels Fig. 2(m), (n), for 
example, the "speed" of the doublon on the light cone 
increases somewhat less than a factor two although J' 
is exactly twice as large. Looking at Eq. (9), this hints 
to a breakdown of the effective model with U — J ^ 0. 
In fact, for U = degenerate perturbation theory in J 
must be considered. Since the states with two fermions 
at the same and at neighboring sites have the same un- 
perturbed energy U = V, decay and recombination of 
the doublon becomes a very efficient process. This leads 
to a maximum mobility as can be seen in panels Fig. 2(i) 
and (o). 

For = V, first-order perturbation theory in J par- 
tially lifts the degeneracy. Therefore, the resulting effec- 
tive model actually describes the motion of a new eigen- 
mode which is a linear combination of a doubly occupied 
site with states where the two fermions are found at adja- 
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cent sites. Rather than doublon propagation, the phys- 
icahy adequate picture is given by propagation of this 
extended object which we wih refer to as an "extended 
doublon" in the fohowing. 

A description by means of an effective model that pre- 
serves the total double occupancy must break down for 
U = 0. This explains the qualitatively different propa- 
gation patterns in the first row of Fig. 2. For V = the 
pattern is given by {Vf^^{t)) = J^_^^{2Jt). For finite V, 
see Fig. 2(e) for example, we note that besides the usual 
propagation pattern describing the delocalization of the 
doublon initially prepared at iq = 50, there is a finite 
probability to find a doubly occupied site around i = 1 
at t ~ 30 J~^. The structure further evolves in time an 
interferes with the main structure. This must be consid- 
ered as a finite-size effect resulting for U = from the 
very fast decay of the doublon into two independently 
moving fermions. Due the periodic boundary condition, 
this implies that the two fermions meet again and form 
a doubly occupied site. The corresponding probability 
strongly decreases with the system size L. 
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V. DECAY OF A DOUBLON AT SHORT TIMES 

The main idea behind the concept of a repulsively 
bound pair of fermions is that energy conservation pre- 
vents the decay of a doublon at strong coupling U > 2W: 
The doublon energy of the order of U cannot be trans- 
ferred to two independently moving fermions with a ki- 
netic energy of the order of at most W each. In Fig. 2, 
the small top panels show the time dependence of the 
total double occupancy. In all cases we find a relaxation 
of the total double occupancy from its initial unit value 
to a nearly constant value after a short time. In many 
cases, this quick initial decay is hardly resolved on the 
scale of the figure, see U = 5 and U = 10 {V = 0) for 
example. 

Fig. 3(a) shows {V{t)) for F = and different U on 
a much shorter time scale up to a few inverse hoppings 
1/J. To quantify the time scale for the doublon decay, 
we look at the position of the first minimum. This "de- 
cay time" is shown in Fig. 3(b) as a function oil/U. We 
find a simple linear relation. The depth of the first min- 
imum also increase with increasing interaction strength. 
In all cases, however, the double occupancy does not re- 
cover completely to its initial value but after some oscilla- 
tions relaxes to a nearly constant value which is becomes 
smaller for weaker U . 

The question how the observed doublon decay is con- 
sistent with energy conservation, is easily answered by 
means of time-dependent perturbation theory in J. For 
J = 0, the total double occupancy is conserved. This 
already explains the high and nearly constant {V(t)) for 
very strong U (see the result for U = 40 in Fig. 3(a)). 
For strong but finite U first-order-in- J time-dependent 
perturbation theory predicts the transition probability 
between two unperturbed energy eigenstates |'0rn) and 



FIG. 3. (a) Short-time behavior of the total double occupancy 
for U = 10, 12, 15, 20, 25, 30, 40 (from bottom to top) and V = 
0. (b) First local minimum point ("decay time") of the time- 
dependent total double occupancy plotted against 1/U. The 
line is a guide to the eye. 

iV^n) to behave as^^ 

9 sin^ (^E^n^A 

i(v^nie-^v^)f a ";/ ^ ;K (11) 

This reminds us that "energy conservation" as used in 
the argument given at the beginning of the section holds 
in the long-time limit only where the r.h.s of Eq. (11) 
evolves into a (5-function. 

Doublon decay is possible (i) at short times or (ii) at 
long times and consistent with energy conservation in 
the presence of additional degrees of freedom to dissi- 
pate the excess energy. Let us discuss the case (i) first 
(see Sec. IX for point (ii)): As a function of the energy 
difference AE^^n^ the transition probability has a peak 
structure with a width that scales as 1/t. Hence, transi- 
tions are possible between states with energy difference 
AEjn^n ^1/^- To put it in other words, excitations 
with energy AEm^n most probably occur on a time scale 
t < 1/AEi^j. Therefore, since the dissociation of two 
fermions in the strong-coupling regime involves energies 
of the order of [/, the position of the first minimum must 
scale with 1//7, as demonstrated in Fig. 3(b). 

At very short times, the decay is independent of the 
coupling [/, as seen in Fig. 3(a) for t < 0.2. This is easily 
explained by Taylor expansion in t: 

{V{t)) = l-t^ AEini + 0{t^) , (12) 

where the variance of the total energy in the initial state 
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is proportional to the number of nearest neighbors z = 2, 

AEini = (V'inil^'l^/'ini) " {'>Pini\U\i'inif = 2zJ^ , (13) 

and thus depends on the hopping amphtude J only. 



VI. DOUBLON DYNAMICS AND 
APPEARANCE-POTENTIAL SPECTROSCOPY 

The time-dependent expectation value of the double 
occupancy at site i is 

mt)) = (0|4e^«*d>;,-«*Jtjo) = \{0\d,e-'»'dl\0)\\ 

(14) 

if a doublon has been prepared at t = at the site zq, 
i.e., dl\0) = cl^cl^\0). Note that the original expectation 
value can be written as a square since (i) H commutes 
with the total particle number and (ii) we start from the 
Fermi vacuum. Namely, starting from the vacuum state 
|0), preparing of the doublon at site zq, time propagation 
and finally annihilation at i, we must return to the same 
state |0). 

The Fermi vacuum corresponds to an empty band in 
the context of electron spectroscopy. Let us discuss the 
relation of doublon dynamics to appearance-potential 
spectroscopy (APS),^^^^ in particular. Consider the fol- 
lowing retarded two-particle (two-electron) Green's func- 
tion: 



For t > we have {Vi{t)) = \Gii^iQi^{t)\'^ and thus 



Gii 



,(t) = -ie(i)(0M; 



|0) . (15) 



This is a ground-state quantity, the Fourier transform 
of which, Giijj{uj + zO+) = ((c-^c-^; c]^c]^))c^, yields the 
appearance-potential spectrum ^^^^^^(a;) = —lmGii^ii{uj-\- 
i0+)/7r.^^ Aii^ii{uj) describes the cross section in a non- 
radiative two-electron process where an initial electron at 
high kinetic energy occupies an empty state in the valence 
band of a metal by transferring the energy difference to 
a core electron which is lifted to another empty state 
in the band. The process is essentially local and repre- 
sents the "time-inverse" of high-resolution CVV Auger- 
electron spectroscopy. 

For an empty band, the equation of motion for the 
APS Green's function is readily solved:^^ 

(16) 
(17) 



with 



uj- e{p) - e{k - p) ' 

Here Ri denotes the position vector to the site z, is a 
wave vector of the first Brillouin zone, and the dispersion 
of the tight-binding band e{k) = —J^^ex.p{—ikA) is 
obtained as a sum over nearest-neighbors displacement 
vectors A. 



(AW) 



27r 



duJe-'^'Gu,^o^o{^^^^^) 



(18) 



At = this related to the Bessel function, {Vf^^{t)) = 
J^ti^{2Jt). For > 0, and using the fact that the 
Green's function is the Hilbert transform of the spectral 
function, we find: 

{Vi{t)) = JdojJ doj'e'^'^-'^'^'Ai^i^,u{co)Au,i,i„{co'). 

(19) 

After substituting uj uo -\- uj' ^ we see that the time 
dependence of the local double occupancy is given by 
the Fourier transform from frequency to time representa- 
tion of the self-convolution of the APS spectral function. 
This relation is remarkable as it provides a link between 
the APS spectral function, an equilibrium quantity de- 
scribing two-particle excitations within the framework of 
linear-response theory, and the non-equilibrium time evo- 
lution of the local double occupancy. It is by no means 
general, however, and can be traced back to Eq. (14) 
which holds in case of an empty band only. 



VII. DECAY OF A DOUBLON 
STABILITY 



LONG-TIME 



In the long-time limit, for an infinitely large system, 
i.e. L ^ 00, the local double occupancy {Vi{t)) for 
t ^ 00 due to a complete delocalization of the doublon 
or the two independent fermions, respectively. The total 
double occupancy (P(t)), however, may relax to a finite 
value. Still there are temporal fluctuations of {V{t)) as 
[D^l-L] 7^ 0, see Fig. 3(b) and also Fig. 2, for examples. 
However, the fluctuations can be quite small as compared 
with the time average 



lim — 

T^oo T 



dt {V(t)) 



(20) 



To quantify these observations, the time average after 
the initial decay at short times as well as the relative 
standard deviation (P^ — V )^/^/P, as a measure for 
the temporal fluctuations, are shown as contour plots in 
Fig. 4. Some sectional views are given in Fig. 5. 

For vanishing couplings U and V the doublon decays 
on a short time scale and is found anywhere in the lattice 
with a probability of approximately 0.019 (for L = 100) 
at later times but fluctuations are strong. 

For finite and increasing /7, but keeping V = 0, the 
doublon stability rapidly rises while the relative fluctu- 
ations decrease. This is understood easily as the energy 
conservation described by Eq. (11), becomes strict in 
the long-time limit, i.e. a single doublon in an otherwise 
empty band is completely stable. 



>^ 




FIG. 4. Long-time average (top) and relative standard de- 
viation (bottom) of the total double occupancy {T>(t)) for 
interaction strengths -10 < [/ < 10 and -10 < F < 10. 
V and (P2 _ V Y^^/V are calculated for the time interval 
50 < t < 100. The color or grey-scale code is given on the 
right. 



Using Eq. (19), we find the time average for T ^ oo, 
^ oc ^ / dujAi^i^^ii{uj)Aii^i^i^{uj) , (21) 

to be given by the integrated square of the non-local APS 
spectral density. As ^^^^^^^^(a;) consists of a finite num- 
ber of J-peaks for any finite L, the integral in Eq. (21) 
is ill-defined. However, one can also compute V directly, 
starting from Eq. (14), inserting resolutions of the unity 
in the form 1 = where |m) is the m-the 

eigenstate of H. Assuming the energy spectrum to be 
non-degenerate and assuming that there is relaxation at 
all, we easily find 



^ = EEKoi4Hl'lmHi^ 



(22) 



With the expressions Eqs. (16) and (18) for the corre- 
sponding Green's function, and using its Lehmann rep- 
resentation, this is also seen to be consistent with Eq. 
(21). Eq. (21) provides the long-time "thermal" value of 
the total double occupancy. 

For large [/, the numerical results of Fig. 4 can be 
perfectly fitted by 



V 1 



m 
IP 



(23) 



with the constant m > as a parameter. This behavior 
can be understood by perturbative arguments using the 





FIG. 5. Sectional views of the stability map Fig. 4 for 
[/, y = 0, 5, 10. Blue lines show P, error bars indicate the ab- 
solute standard deviation. Green lines with error bars: time 
average and standard deviation of the nearest-neighbor oc- 
cupancy, (I E<ij> Ectct' ^*v(^)^^cT'(^))- Red: sum of double 
and nearest-neighbor occupancy. 



canonical transformation discussed in Sec. III. Using Eq. 
(6), with the generator 5, 



-iS A'H.' t AS 



-iS —il-L' t AS j] 



e^^dlM, (24) 



and exploiting particle- number conservation, we get: 

m)) = Y.\m,e-'^e'^''e'^d\\Q)'\ . (25) 

i 

The state 

m ^ e'^4\0) = dl\0) + iiSJlm + 0{jyu') (26) 

is a linear superposition of a one-doublon and a zero- 
doublon state: 

K) = K,i)-^K,o) + ^(^V^'), (27) 



with li^l,) = dl\0) and |V^,o) = 

cj^c^j^)|0). This characteristic is, separately, preserved 

under the time evolution e'^ *. After some algebra we 
find: 



h.c. + 0( jVf/4) 



(28) 



Hence, perturbative in 1/U corrections to the total dou- 
ble occupancy are of the order J'^ /U'^ . 
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Eq. (28) furthermore shows that in leading perturba- 
tion order there is a separation of characteristic time 
scales. The first matrix element involves energies in the 
one-doublon subspace and thus a short time scale 1/U. 
This is the time scale of the strong initial oscillations of 
{V(t))^ seen in Fig. 3, and once more explains the 1/U de- 
pendence of the "decay time" , i.e. the position of the first 
minimum. The second matrix element between states in 
the zero-doublon subspace provides a longer time scale 
~ 1/ J. This is the scale on which the oscillations decay. 
Finally, corrections to the 1/U scale are provided by ef- 
fective doublon-hopping processes. This results in a scale 
1/ J' ~ the effects of which, however, are too weak 

to be seen in Fig. 3. 

Finally, let us discuss the results for a finite nearest- 
neighbor interaction V . As is shown in Figs. 4 and 
5, the initially prepared doublon is most unstable for 
V ^ U. Here, the transition d\ |0) 4t4±i,i 1^) 
comes resonant. A further separation of the fermions 
beyond nearest-neighbor distances, however, is the more 
suppressed the larger V gets. The latter is obvious for 
reasons analogous to those given above regarding the U 
dependencies. As already noted in the context of prop- 
agation patterns above, ioi U = V ^ 0^ an "extended 
doublon" is formed as a linear combination of a dou- 
bly occupied site with states where the two fermions are 
found at adjacent sites. Though the probability for find- 
ing two fermions at the same site anywhere in the lattice 
shows a minimum for U = V in the stability map in 
Fig. 4, the one for finding them as nearest neighbors is 
almost equally large as can be seen in the sectional views 
of the stability map in Fig. 5. Furthermore, the sum of 
both equals the one for the same value of U but vanishing 
V. At the same time the oscillations of the double occu- 
pancy as well as the one of nearest neighbor state exhibit 
a maximum (Fig. 4 and 5) whereas their sum does not. 
Hence the oscillations cancel each other. 



VIII. TWO DOUBLONS 

The preceding examinations were restricted to the sub- 
space of two fermions. This lacks some important as- 
pects, such as doublon-doublon and doublon-fermion- 
scattering. In the following we therefore extend our study 
to four-fermion states. The size of the one-dimensional 
lattice is fixed to L ^ 50. 

To begin with, consider an initial state at t = with 
two doublons at neighboring sites: |V^ini) = dldj |0) with 
\i — j\ = 1. In the strong-coupling limit U^V ^ J, 
this state has a mean energy of the order of 2/7 + 4V + 
0{J'^ /U^ J'^ /V): A state with two neighboring doublons 
entails two neighboring fermions for each constituent 
fermion. Processes starting from this state and involving 
a single or two hopping events will dominate the physics 
in the strong-coupling case and are sketched in Fig. 6. 

Fig. 7 shows the time-dependent local and total double 
occupancy for different U and V. The overall trends can 




FIG. 6. Scheme of dominant first-order [(1), right] and 
second-order [(2), left] hopping processes from an initial state 
with two doublons on nearest-neighbor sites to the possible 
final states. The involved interactions U and V are depicted 
by wiggly lines in red and blue, respectively. 



easily by understood by focussing on certain resonant 
cases: 

(i) If V = —U the first order process, referred to as (1) 
in Fig. 6, becomes resonant: The initial and final state 
have the same mean energy up to a small correction of 
the order 0( J^/[/, J^/F). In the strong-coupling limit a 
further spatial separation of the fermions is suppressed 
as there is a large excess energy U or V that cannot be 
accommodated in the system. A propagation of the com- 
pound object over many lattice sites is only possible via 
second-order hopping processes with a very low probabil- 
ity as compared to the first-order process (1). We there- 
fore expect the two doublons to be basically localized at 
their initial positions. This explains the pattern shown 
in Fig. 7(ann). 

After some settling time the total double occupancy 
(see top panel in Fig. 7(ann)) tends to a value slightly 
less than unity which is less than expected for both states 
that define the process (1). We therefore conclude that 
there is a certain non-zero probability for the decay of 
the compound object into fragments without double oc- 
cupancy that is not consistent with energy conservation. 
As discussed for the two-fermion case, this is possible at 
very short times. 

The main dynamical effect, however, consists in a rapid 
oscillation between the two states of process (1). In the 
map for the time average P, see Fig. 8 (left), this man- 
ifests itself as a "valley" along the bisecting line of the 
second and fourth quadrant. Furthermore, this is accom- 
panied by a maximum in the relative fluctuations (not 
shown), similar as in the two-fermion case. 

(ii) Correspondingly, we find another "valley" along 
the line given by V = —2/7 in Fig. 8 (left). This is 
associated with the second-order process (2^) in Fig. 6 
which becomes resonant in this case. Again, there is 
mainly an oscillation between the two states of (25) which 
both have the energy 2/7 +4V = 3^ . The process involves 
a virtual intermediate state with an off-resonant energy 
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0.55 



(asep) U = 5,V = 5 




(Cnn) U = 0,V = 10 



(dnn) U = 5,V = 



(Cnnn) U = 0,V = 10 



(dnnn) U = 6,V ■ 



(bsep) U = 10,V -- 




(enn) f/= 10,y : 



(fnn) f/ = 10,y = 10 



(ennn) U 



(fnnn) f/ ^ 



(Csep) U = 10,V = 



FIG. 7. Time dependence of the expectation value of the local (main panels) and the total double occupancy (small top 
panels) for different U and V as indicated. Axis and color code as in Fig. 1. Calculations are performed for different initial 
states as indicated by the bracketed symbols: two pairs fermions (two doublons) initially prepared as nearest-neighbors (xnn), 
next-nearest-neighbors (xnnn) or further separated with = 10 (xgep) respectively. Results for one-dimensional lattice with 

periodic boundary conditions and L = 50, 51 and 49 sites, respectively. 



U^3V. 

As before in case (i), a propagation of the compound 
object over many lattice sites is suppressed as it necessar- 
ily involves fourth-order processes. In fact, Fig. 7(bnn) 
shows that the fermions essentially remain close to their 
initial sites. 

An oscillation between the two states of (2^) clearly 
implies the total double occupancy to oscillate between 
approximately 2 and 0. In the long-time limit it tends to 
relax to a value close to or slightly less than 1. 

(iii) In case of vanishing [/, the second-order process 
{2a) becomes resonant at the energy 2/7 + = 4V. 
This causes another branch of minima along the F-axis 
in Fig. 8 (left). 

Opposed to cases (i) and (ii), the four-fermion cluster 
may propagate via the (2^) process followed by a pro- 
cess inverse to (2^) but resulting in two neighboring dou- 
blons shifted by one site to the left or right as compared 
with the initial state. Repeated second-order hopping 
processes then lead to a more efficient delocalization of 
the cluster and thus also of the expectation value for the 
double occupancy as is seen in Fig. 7(cnn)- 

(iv) For a vanishing V, the process (2^^) becomes res- 
onant at the energy 2U + — 2U. This implies 
that the initial cluster with two neighboring doublons 
can dissociate into two doublons separated at arbitrar- 



ily large distances via second-order hopping processes 
over off-resonant intermediate states. Delocalization is 
thus very efficient and results in the pattern displayed in 

Fig. 7(dnn). 

The propagation pattern is obviously dominated by 
two "light cones" with different velocities. This can 
be traced back to the interaction between the two dou- 
blons by comparing with the patterns in Fig. 7(bsep) and 
Fig. 7(csep) which refer to an initial state where the two 
doublons are well separated and prepared at a distance 
= 10 and where the mode with lower velocity is ab- 
sent. It is an open question whether the slow mode is due 
to the repulsive hard-core constraint or due to the attrac- 
tive interaction in the effective Hamiltonian Eq. (8). The 
"light cone" associated with the higher velocity is identi- 
cal to the one found for propagation of a single doublon, 
see Fig. 7(dnn) and Fig. 2(h) and mind the different lat- 
tice sizes. 

In Fig. 8 (left), we find a signature of the resonant 
process (2^^) along the = line. As in the two-fermion 
case, the doublons are stabilized with increasing Hubbard 
U. 

(v) Finally, the process (2c) gets resonant if 2U -\-4y = 
U -\-2V which again becomes manifest in a "valley", given 
by y = — 1/7, in the map. Fig. 8 (left), which is clearly 
visible at larger values of U and V. 
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(a) nearest-neighbors (b) next-nearest-neighbors (c) further separated 

FIG. 8. Time average of the total double occupancy {T>{t)) for interaction strengths — 10 < [/ < 10 and — 10 < V < 10. 
Calculations for an initial state with two doublons placed at sites i and j. (a) i and j nearest neighbors, (b) next- nearest 
neighbors and (c) |i — i| = 10. Average over the time interval 50 < t < 100. The color code is given on the right. 



Regarding the mobility, we note that the process (2c) 
can be either inverted or the fermion triple can move 
resonantly through the lattice. Both possibilities con- 
tribute to the computed propagation pattern shown in 
Fig. 7(enn). 

In all other cases, the initial state shows both a high 
stability and a marginal mobility in the strong-coupling 
limit. Fig. 7(fnn) gives an example for /7 = V = 10. 
We note that the relative fluctuations around the time 
average amounts to approximately 1% only. 

Although the underlying physics is the same, the re- 
sults are completely different if the two doublon are pre- 
pared at sites which are next-nearest neighbors. The 
calculated propagation patterns are shown Fig. 7 in the 
third and fourth column, while Fig. 8 (middle) displays 
the corresponding time averages. The figures also show 
the results for two doublons prepared initially at a dis- 
tance \i — j\ = 10 (Fig. 7, last column and Fig. 8, right). 

For two doublons prepared as next-nearest neighbors, 
the dominant first-order and second-order hopping pro- 
cesses are sketched in Fig. 9. First, we note that the 
processes (Ij,), (2^), (2^) and (2j) are all independent 
of the problem's four-particle character. Provided that 
the physics is dominated by those processes, one would 
expect the propagation pattern of two initially next- 
nearest-neighboring doublons to essentially resemble that 
of two independent doublons. In the strong-coupling 
limit, this is the case for processes (Ij,), (2^) if U = V 
and independently of V for (2j). As is seen Fig. 7(annn) 
and (bnnn), the doublons' propagation is described by 
about the same maximal effective hopping as in the case 
of a single doublon, see Fig. 2(o), for example. There is, 
however, an additional mode visible in Fig. 7(annn) and 
(bnnn) which rcsults from the two doublons resting more 
or less at their initial sites. This is caused by the respec- 
tive inverse hopping processes and basically disappears 
with increasing interaction strengths U = V ^ oo and 
also in the case where the two doublons are prepared at a 
larger distance (see Fig. 7(asep)). A branch of minima oc- 
curs along the line U = V in the stability map, Fig. 8(b), 
which looks similar to that obtained in the two-fermion 



case (cf. Fig. 4). The process (2^) is resonant only if 
U = 0. Here the doublons rapidly dissociate into more 
or less independent fermions resulting in deep "valley" 
around U = in Fig. 8(b). 

The processes (1[^), (2^^), (2J,) and (2^^) are immanent 
to the four-particle character of the problem and become 
resonant if V = ^U, V = \U , V = \U or V = ^U, 
respectively. The same holds for the inverse to process 
{2d) (see Fig. 6) which becomes resonant if V vanishes. 
Except for the last one, the doublon number is changed 
in all processes. We therefore expect and find a region 
of instability, bounded from below by V = as can 
be seen from the level curves in Fig. 8(b). Generally, 
the propagation patterns 7(dnnn)7 (ennn) and (fnnn) are 
not easily interpreted by means of simple perturbative 




FIG. 9. Scheme of dominant first-order [(1^5 right] and 
second-order [(2^), left] hopping processes from an initial state 
with two doublons on next-nearest-neighbor sites to the pos- 
sible final states. The inverse of process (2d) (see Fig. 6) is 
not shown again. U and V are depicted by wiggly lines in red 
and blue, respectively. 
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arguments. 



IX. DOUBLON-FERMION SCATTERING 



It is worth mentioning that for vanishing nearest- 
neighbor-interaction V = (not displayed) the doublons 
essentially show the same spreading behavior as they did 
in the case of a single doublon (see Fig. 7(dnn)) and their 
stability again rises with \U\. Further, for large cou- 
plings of opposite sign U = — V, all processes except 
for (2j) are strongly suppressed. The patterns (not dis- 
played) are rather similar to those for a single doublon 
(see Fig. 7(ann)). 

The further away two doublons are prepared in the ini- 
tial state the less they influence each other. We therefore 
obtain results similar to those for a single doublon. This 
can be seen from our calculations with two doublons ini- 
tially separated by 10 sites by comparing e.g. the maps 
for the long-time averages P, Fig. 8(c) and Fig. 4, as well 
as by comparing the propagation patterns in Fig. 7 and 
2 for corresponding interaction strengths. 

Generally, the propagation patterns considerably dif- 
fer from the corresponding ones for doublons formed by 
bosons. Motivated by experiment,^ Petrosyan et al.^'^ 
consider the Bose-Hubbard model, 1-L = — ^Xl<ij> + 
{U/2) n^{n^ — 1)7 in the strong-coupling limit with an 
additional constraint excluding states, analogous to the 
Fermi case, with two or more bosons at the same site. 
Preparing an initial state with two neighboring doublons, 
propagation patterns are obtained which look very sim- 
ilar to our cases U = —V = 10 or U = V = 10 (see 
Figs. 7(ann) and 7(fnn)), i-e. propagation is strongly sup- 
pressed. This can be understood by again referring to a 
respective effective model for the strong-coupling limit. 
Canonical transformation yields:^'^ 



K 



eff ^ 



<ij> 



Here, d^^^ denotes the annihilation (creation) operator ent initial states. Besides d]^ 11^2), we also study the 



2 J' ■ 



<ij> 



(29) 



The propagation and the decay of a repulsively bound 
pair is expected to be strongly affected by the presence of 
additional fermions. As a finite fermion density cannot be 
studied reasonably by means of the Krylov approach, we 
will here consider two additional fermions only. To this 
end we first determine the ground state of the Hamilto- 
nian in the two- fermion subspace 11^2) and subsequently 
add a doublon at a certain site io to define the initial 
state Since the weight of doubly occupied sites 

in the ground state is almost vanishing for a lattice with 
L = 50 sites, this setup allows to study the scattering 
of the doublon with almost independently propagating 
fermions. 

Fig. 10 demonstrates that there is in fact a strong im- 
pact of the additional fermions on the doublon prop- 
agation. This crucially depends on V. For V = 
the propagation pattern is similar to the one shown in 
Fig. 2(m). For finite however, we notice the pres- 
ence of essentially two modes. This is most obvious for 
U = V. Comparing with the results in the absence of ad- 
ditional fermions, i.e. Figs. 2(n) and 2(o), the fast mode 
can be interpreted without referring to the presence of 
the fermions as before: For U = 2V there is a reso- 
nant second-order hopping process moving the doublon 
one site further. For U = V even a first-order process 
is resonant which explains the higher mobility in that 
case. There is, however, an additional slowly propagat- 
ing mode which in the case U = V even carries almost all 
weight. As this mode is absent in Figs. 2(n) and 2(o), it 
is attributed to fermion-assisted doublon propagation or 
rather a suppression of the mentioned resonant processes 
due to the presence of additional fermions which can be 
discussed in a way similar to the two-doublon case above. 
One should note, however, that the slow mode must dis- 
appear for larger lattices as the probability for doublon- 
fermion scattering events tends to zero for L ^ 00. 

As concerns the decay of the doublon we therefore fo- 
cus on the V = case only and also consider differ- 

■f 



for doublons made up of bosons {U^"^). As in the Fermi 
case, the effective hopping is given by J' = AJ'^ /U . Eq. 
(29) should be compared with Eq. (8). In contrast to 
the fermionic case, the attractive interaction between two 
nearest-neighboring doublons is larger by a factor four for 
doublons made of bosons. This explains the tendency to 
a strongly suppressed propagation. 

It also explains that, in the bosonic case, the formation 
of clusters of doublons is favored and phase separation 
is possible below some critical temperature.^'^ Contrary, 
in the Fermi case, doubly occupied sites may Bose con- 
densate under certain circumstances.^^ In fact, we did 
not find any indications for a clustering of doublons. 
Two doublons are rather never found to form a bound 
state unless an explicit nearest-neighbor-interaction V is 
present. 



system's time evolution starting from states where two 
fermions are prepared at sites close to the initial posi- 
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(a) [/ = io,y : 



(b) [/ = io,y = 10 



FIG. 10. Time evolution of the local and total double oc- 
cupancy when a doublon was initially created in the ground 
state of the two-particle system. The presentation is the same 
as in Fig. 1. 
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FIG. 11. Coefficient m obtained from a fit of the time- 
averaged expectation value for the tot al do uble occupancy 
V{t) to the observed U dependence: V(t) (^^(0))(1 - 
m/U'^). m values are obtained from linear regression of our 
data in the range 8 < t/ < 45. Results are shown for different 
initial states as indicated and discussed in the text, n = 1 
refers to the results of a time-dependent DMRG calculation 
by Al-Hassanieh et al., see Ref. 18 and text for discussion. 



tion of the doublon zq, i.e. |m, m') = 4o+mt<l4o-m' JO)- 
This is compared to results obtained for two doublons at 
nearest-neighboring sites, <^1q<^1q+i |0), and two doublons 
prepared at a distance of 2 and 10, i.e. <^Jq<^Jq+2|0) 
^^Jq^^Iq+ioIO)? respectively. In all cases we find a decay of 
the doublon expectation value on a short time scale 1/U 
followed by a stabilization to a nearly constant value at 
large times. The residual quantum fluctuations are dis- 
regarded by looking at the time average V. As before, 
we find that the decayed doublon fraction scales linearly 
with l/[/2 for large times, V(t) {V{0)){1 - m/U^). 
Hence, in the strong-coupling limit the doublon stability 
is quantified by the coefficient m. For m = there is no 
decay at all, and a small value for m indicates a rather 
stable doublon. Our results for the different initial states 
are shown in Fig. 11. 

Generally, for a system with additional fermions, one 
expects an hugely increased phase space for inelastic pro- 
cesses leading to doublon decay. On the other hand, the 
energy-conservation argument suggests that for strong U 
a rather complex inelastic process has to take place to 
allow for decay, namely a process of high order where a 
sufficient number of particles must be involved to dissi- 
pate a large energy of the order of U. While such pro- 
cesses are expected to be exponentially suppressed for 
large [/, they should contribute to some degree and lead 
to a destabilization of a doublon. 

However, our results for different initial states, as dis- 
played in Fig. 11, just show the opposite trend: The 
presence of two additional fermions in the initial state in 
all cases leads a smaller coefficient m in the decay 
law. The strongest effect is visible for the initial state 



|1, 1) where the two fermions are neighbors of the dou- 
blon at iQ. Here m is the smallest and the doublon is 
most stable, m increases with increasing distance of one 
of the fermions from the position of the doublon, see the 
initial states |1,2) and |1,3). It further increases if also 
the second fermion is positioned at a distance from io 
(see |2,2), |3,3) and |4, 4)), and it approaches the value 
obtained for the case where both fermions are delocalized 
in the ground state d]^ 1^2)- The maximum value is ob- 
tained for the isolated doublon in an otherwise empty lat- 
tice, i.e. for dljO). If the two fermions themselves form a 

second doublon, see the results for dl^dl^_^JO) in Fig. 11, 
this again tends to stabilize the original one: m decreases 
with decreasing distance x between the two doublons. 

These trends can be understood if the doublon dynam- 
ics is considered at short times: First-order-in- J time- 
dependent perturbation theory shows that doublon decay 
is allowed on a time scale 1/U as has been detailed in Sec. 
V. Here, one can argue that an unoccupied site neighbor- 
ing the doublon is necessary for the decay process as the 
immediate surrounding is relevant for its start. Hence, 
the coefficient m is the smaller and the doublon is more 
stable if decay channels are blocked by localized fermions 
or doublons close to the doublon at io and, to a lesser 
extent and depending on the size of the lattice, even by 
two delocalized fermions in the two-fermion ground state. 
This nicely explains the results described above. 

After that time scale, energy conservation as expressed 
by Fermi's golden rule, applies and the total double occu- 
pancy virtually relaxes to a constant value. As analyzed 
in Sec. VII, the probability for the dissociation of a dou- 
blon should then scale as 1/U'^. On an for large U ex- 
tremely long time scale, which exponentially depends on 
[/,^^ contributions from higher-order perturbation theory 
in J/U become important and would generally allow for 
further decay in more complex processes. 

In this context it is interesting to compare our re- 
sults with the those of a time-dependent density-matrix 
renormalization-group (DMRG) study by Al-Hassanieh 
et al.^^ where the decay of a doublon created by a 
nearest-neighbor particle-hole excitation of a half-filled 
one-dimensional Fermi Hubbard model was considered. 
The DMRG calculations show (i) a fast decay at a char- 
acteristic time scale l/[/, (ii) a basically constant double 
occupancy at larger times up about 40 J~^, and (iii) a 
scaling of the decayed fraction of the doublon. All 
this agrees perfectly with our results obtained for four 
fermions only. The m coefficient taken from the DMRG 
results^^ is also included in Fig. 11 ("n = 1") and is 
found to be close to that obtained for the |1,2) initial 
state. Even this is plausible since the spin-dependent 
site occupations of the state |1, 2) and of the initial state 
of the DMRG calculation are the same in the immediate 
environment of io. Note that a quantitative compari- 
son with the DMRG study of Ref. 21 for the half-filled 
Hubbard model is not possible as an initial state where 
doubly occupied and empty sites alternate, is considered 
there. Still the qualitative features are rather similar. 
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We conclude that there is at least a qualitative agree- 
ment with the dynamics of the half-filled model on the 
time scale accessible to time-dependent DMRG. Hence, 
the local blocking of decay channels determines the dou- 
blon decay on this time scale while an exponential decay 
of the doublon, associated with the imaginary part of the 
doublon-propagator at frequencies ~ [/, may take place 
on a very much larger time scale only. At short times, the 
stability of a doublon is thus enhanced by the presence 
of additional fermions while a destabilization, due to an 
increased phase space for decay, appears irrelevant. 



X. SUMMARY 

Concluding, the real-time dynamics of two or a few 
more strongly interacting Fermions moving in a peri- 
odic lattice potential exhibits a surprisingly rich physics 
which is not only linked to experiments with ultracold 
atoms trapped in optical lattices but also to electron 
spectroscopy of metal surfaces as well as to rather gen- 
eral questions on the propagation and decay of bound 
quantum states and the relaxation of quantum systems 
prepared in a highly excited initial state. Unfortunately, 
there is a rather limited number of theoretical tools to 
analyze quantum dynamics with a macroscopically large 
number of fermions far from equilibrium. Either one has 
to tolerate mean-field type approximations like in the 
non-equilibrium dynamical mean-field approach^^^^ or 
has to restrict oneself to one-dimensional or impurity- 
type systems to render an application of time-dependent 
renormalization-group approaches^^^^ possible. Exact- 
diagonalization or Krylov-space methods^^^^ are numer- 
ically exact and fiexible but strongly restricted to few- 
particle systems with a moderately large Hilbert-space 
dimension only. 

Even the analysis of the two-fermion case, however, 
helps to understand important concepts such as the tem- 
poral stability of a doublon, i.e. a repulsively bound 
pair of fermions. The decay of a doublon in an oth- 
erwise empty system is possible on a very short time 
scale given hy 1/U where energy conservation, within the 
spirit of time-dependent first-order perturbation theory 
or Fermi's golden rule, does not apply. Using pertur- 
bative diagonalization of the Hamiltonian by means of 
a canonical transformation, one can understand the ob- 
served dependence of the fraction of the doublon 
that has decayed in the long-time limit. The time aver- 
age of the total double occupancy is found to be given 
by a quantity defined for the equilibrium or ground state 
of the system, namely the integrated square of the spec- 
tral density related to appearance-potential spectroscopy. 
But also the fully time-dependent local double occupancy 
can be expressed in terms of this spectral function, which 
must be seen as an unexpected interrelation valid for a 
two-particle system only. 

The spatiotemporal evolution of the expectation value 
of the local double occupancy can be expressed in terms 



of Bessel functions of the first kind for = as well as 
in the strong-coupling limit. For the case of a non-zero 
nearest-neighbor interaction U, the propagation pattern 
is a bit more complicated but can be understood with 
the help of perturbative arguments again. We find that 
depending on V the propagation mechanism may change 
qualitatively. Upon approaching the case U = the 
propagation becomes more efficient as a first-order-in- J 
rather than a second-order process becomes resonant. 

In the case of four fermions, the propagation pat- 
terns are much more complicated. Referring to first- and 
second-order hopping processes which, depending on the 
values of the interaction strengths U and U, do or do 
not become resonant, we have analyzed the real-time dy- 
namics after preparation of different initial states, such 
as two neighboring doublons or doublons prepared at a 
certain distance from each other as well as a single dou- 
blon on top of the two-fermion ground state or with two 
fermions localized at distinct distances from the doublon 
in the initial state. Even for the strong-coupling limit, 
we found that the temporal evolution of the local dou- 
ble occupancy can be understood in most but not in all 
cases by perturbative arguments. The physics of a fi- 
nite density of doublons consisting of fermions is known 
to be rather different from the case of doublons made of 
bosonic particles which undergo a transition to a phase- 
separated state instead of Bose condensation.^'^' ^'^ Con- 
sistent with this, we did not find any indications for a 
clustering of doublons consisting of fermions unless an 
explicit nearest-neighbor-interaction V is present. 

Surprisingly, there is a rather regular trend concerning 
the decay of a single doublon in the presence of two more 
fermions. We found that the total double occupancy, 
apart from quantum fiuctuations, relaxes to a constant 
value after an initial decay on a time scale 1/U. The long- 
time average deviates from the initial value by a fraction 
that scales with U as l//7^ in the strong-coupling limit, 
like in the case where there are no additional fermions, 
but with a coefficient m that characteristically depends 
on the initial state, m is found to decrease and thus the 
stability of the doublon is found to increase when two 
fermions are added - a result which at first sight is con- 
fiicting with the expectation that adding more degrees of 
freedom to the system should strongly increase the phase 
space available for decay in energy-conserving processes 
where the doublon energy U is dissipated to a large num- 
ber of particle-hole or spin excitations. Those processes, 
however, require a huge time scale to contribute signifi- 
cantly to the doublon decay. More important for the sta- 
ble fraction of the doublon is the local environment in the 
initial state as the main effect of an additional doublon 
or of additional fermions in its vicinity is to block decay 
channels on the short time scale on which decay is pos- 
sible rather than ruled out by energy conservation. This 
is a general argument which apparently also applies to 
the half-filled case, for example. In fact, we find almost 
quantitative agreement with a time-dependent DMRG 
calculation.^^ On the other hand, the argument leaves 
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the possibility for an e.g. exponential-in-t decay law on 
much larger time scales which might be expected on gen- 
eral grounds. ^^'^^ 
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